Map of Corfu

Library loading

library(tidyverse)
library(sf)
library(ggtext)
library(sp)
library(raster)
library(patchwork)
library(extrafont)
library(osmdata)
library(rosm)

I had to cancel my sailing trip to Corfu - Minor annoyance, minor in the context of a global pandemic disrupting many lives. That is the reason why I chose to augment part of the sailing map I would be using with some additional features such as lighthouses and relief shading.

This entire exercise would not only allow me to project features on maps, but also allow me to project myself into this map… The idea of mixing a “real” picture with relief data came a tweet of Tyler Morgan Wall showing an old map of India.

https://twitter.com/tylermorganwall/status/1281950095756402689?lang=en

I started by taking a picture:

In order to use it, I had to add some georeference data to the picture I followed the following tutotrial using the QGIS application:

https://www.qgistutorials.com/en/docs/georeferencing_basics.html

corfu_map = raster::stack("data/corfu_wood.tif")

# we extract the bounding value from the georeferenced file

 bbx1<-extract_bbox(corfu_map)

#we reduce it's size as the shapefile and elevation data includes
#some part of mainland Greece
min_lon<-19.620209
max_lon<-19.967651

min_lat<-39.6248
max_lat<-39.832796

min_lon->bbx1[1,1]
max_lon->bbx1[1,2]

min_lat->bbx1[2,1]
max_lat->bbx1[2,2]

# we create the bbx tibble and polygon which we will use for cropping

bbx_tb1<-tibble(
  x=c(min_lon,min_lon,max_lon,max_lon),
  y=c(max_lat,min_lat,min_lat,max_lat)
)


bbx_polygon1<-bbx_tb1 %>%
  st_as_sf(coords = c("x", "y"), crs = 4326)%>%
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("POLYGON")


#we load the shapefile and transform it to the desired projection
# and crop it

corfu_shp<-st_read("data/periphereies/periphereies.shp")

## Reading layer `periphereies' from data source `/Users/andrewmacpherson/Github/public/Cartography/cartography/data/periphereies/periphereies.shp' using driver `ESRI Shapefile'
## Simple feature collection with 13 features and 1 field
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 104011.4 ymin: 3850811 xmax: 1007936 ymax: 4623914
## projected CRS:  GGRS87 / Greek Grid

corfu_shp_wsg<- corfu_shp %>%
  st_transform(crs=4326)

cropped_corfu<-st_crop(corfu_shp_wsg,bbx_polygon1)

# we do the same for the raster file
srtm_corfu <- raster("data/corfu_hill.tif")
cropped_srtm_corfu<-mask(srtm_corfu,bbx_polygon1)
cropped_srtm_corfu_spdf <- as(cropped_srtm_corfu, "SpatialPixelsDataFrame")
cropped_srtm_corfu_tb <- as_tibble(cropped_srtm_corfu_spdf)


  corfu_map_2 <- corfu_map%>%
    # hide relief outside of Corfu by masking with country borders
   mask(cropped_corfu, inverse=TRUE) %>%
    as("SpatialPixelsDataFrame") %>%
    as.data.frame() %>%
    mutate(hex=rgb(corfu_wood.1,corfu_wood.2,corfu_wood.3, maxColorValue = 255))


#we extract some features
  highways <- bbx1 %>%
    opq()%>%
    add_osm_feature(key = "highway",
                    value=c("motorway", "trunk",
                            "primary","secondary",
                            "tertiary","motorway_link",
                            "trunk_link","primary_link",
                            "secondary_link",
                            "tertiary_link")) %>%
    osmdata_sf() %>% 
    trim_osmdata(bbx_polygon1) 

  light<- bbx1 %>%
    opq() %>%
    add_osm_feature(key = "man_made" ,
                    value=c("lighthouse")) %>%
    osmdata_sf()%>% 
    trim_osmdata(bbx_polygon1) 

 corfu<- ggplot() +

    geom_raster(data = cropped_srtm_corfu_tb, aes(x = x, y = y, fill = corfu_hill))+
    scale_fill_gradient(low = "black", high = "white")+
    geom_raster(
      data = corfu_map_2,
      mapping = aes(
        x = x,
        y = y), fill=corfu_map_2$hex)+

    geom_sf(data=highways$osm_lines,
            colour="darkred",
            size = 0.3,
            alpha = 1)+
    geom_sf(data=light$osm_points,
            fill="red",
            colour="red",
            size = 1,
            alpha = 1)+
    theme_void()+
 
   labs(title="Corfu, Greece",
        subtitle="Augmented sailing map, elevation, roads and lighthouses",
        caption="source: Imray North Ionian Islands G11 map, https://dwtkns.com/srtm30m/, https://earthdata.nasa.gov, OpenStreetMap")+
   
   theme(
     plot.title = element_text(color = "#010057", size = 25, face = "bold"),
     plot.subtitle = element_text(color = "#010057", size = 15),
     plot.caption = element_text(color = "#010057", face = "italic"),
     legend.position = ""
 )


 corfu